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t^^ . In this paper, three multiscale methods for couphng of mesoscopic (compartment-based) 

^ \ and microscopic (molecular-based) stochastic reaction-diffusion simulations are investigated. 

Two of the three methods that will be discussed in detail have been previously reported 
w ' in the literature; the two-regime method (TRM) and the compartment-placement method 

(CPM). The third method that is introduced and analysed in this paper is the ghost cell 

method (GCM). Presented is a comparison of sources of error. The convergent properties 
^ \ of this error are studied as the time step At (for updating the molecular-based part of the 

^ O^ l model) approaches zero. It is found that the error behaviour depends on another fundamental 

computational parameter h, the compartment size in the mesoscopic part of the model. Two 
K^ important limiting cases, which appear in applications, are considered: 

(^ . (i) At — )■ and h is fixed; 

0\ ; (ii) At -^ and /i -^ such that A/At//i is fixed. 

[^^ \ The error for previously developed approaches (the TRM and CPM) converges to zero only 

— u ' in the limiting case (ii), but not in case (i). It is shown that the error of the GCM converges 

O ■ in the limiting case (i). Thus the GCM is superior to previous coupling techniques if the 

mesoscopic description is much coarser than the microscopic part of the model. 
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1. Introduction 

Multiscale stochastic reaction-diffusion methods which use models with different levels of 
detail in different parts of the computational domain are applicable to a number of biological 
systems, including modelling of intracellular calcium dynamics [l2|, MAPK pathway [20] 
and actin dynamics j9|] . In these applications, a detailed modelling approach (which requires 
simulation of trajectories and reactive collisions of individual biomolecules) is only needed 
in a small part of the computational domain. The main idea of multiscale methods is then 



simple to formulate [10[: we use a detailed modelling approach in the small subdomain 



of interest and a coarser model in the rest of the computational domain. In this paper. 
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detailed molecular-based (microscopic) models will be given in terms of Brownian dynamics 
[sl, |25| . The remainder of the computational domain will be divided into compartments and a 
mesoscopic (compartment-based) model will be used, i.e. we will simulate the time evolution 
of the numbers of molecules in the corresponding compartments p, [l9|. 

There have been a number of approaches developed for coupling different reaction- 
diffusion models. They include coupling of meso scop ic (compartment-based) models with 
coarser (mean- field) PDE-based descriptions [l3l. I2I. l26l. l23l| . coupling of microcopic (molecular- 



based) models with mean- field PDEs |l5l.ll8l.ll4l|. and coupling of microscopic and mesoscopic 
models [lO|, llll, l20|, Isjj A successful multiscale algorithm requires an accurate implementation 
of inter-regime transfer of molecules. In this paper, we will study convergence properties of 
two algorithms for coupling microscopic and mesoscopic descriptions which were previously 



published in the literature: the two-regime method (TRM) [id, lUl] and the compartment- 
placement method (CPM) J20l]. One of the conclusions of our analysis is that these algorithms 
do not converge in the limit of small time steps and a fixed compartment size. Thus, we also 
propose another approach, the ghost cell method (GCM) which is suitable for this parameter 
regime. 

We will consider a react ion- diffusion model in the computational domain Q C M^ for 
both A^ = 1 and A^ = 3. We will divide Q into two parts, open sets Qm and Qc^ which 
satisfy 

n^un^ = n and QA:nQc = ^, (1) 

where an overbar denotes the closure of the corresponding set. The microscopic simulation 
technique is used in Qm- Each molecule, j, in Qm is considered to be a point particle at 
some location in space, Xj(t), at time t, which is updated according to discretized Brownian 
motion, i.e. 



Xj(t + At) = Xj(t) + ^2DjAtC, (2) 

where Dj is the diffusion constant of the j-th molecule. At is a small prescribed time step 
and C is a vector containing zero mean, unit variance normally distributed random numbers. 
In this paper, we will study the convergence of multiscale methods in the limit At — )■ 0. 
Since the discretized Brownian motion ([2]) is only used in Qm, we have to specify what will 
be done in the remainder of the domain, fi^, where the mesoscopic model is used. In this 
paper, we distinguish the following two cases: 

(i) the mesoscopic model is kept fixed in the limit At — ^ 0; 
(ii) the mesoscopic model is refined as At approaches zero. 

The resolution of the mesoscopic model (compartment size) will be denoted by h. Of par- 
ticular interest is the error that is caused as a direct result of the coupling and thus we will 
use the parameter h as a measure of the compartment size at/on the interface between the 
two modelling subdomains. In the case of regular cubic compartments of volume h^, the 
parameter h is simply the length of an edge of each cube. We will also consider unstructured 
meshes where the compartment size h will be suitably generalized. Using h, the cases (i)-(ii) 
can be formulated as follows: 

(i) At — ;■ and h is fixed; 



(ii) At — )■ and h -^ such that yAt/h is fixed. 
Both hinits (i) and (ii) are important in apphcations. We will see that the error at the 



interface d^M H dilc of previously developed methods [lO, lUj, l20|] converges to zero in the 
limit (ii). This limit requires the refinement of the mesoscopic model. However, the standard 
mesoscopic model converges in the limit /i — )■ only if the molecules are subject to zero- 
order or first-order chemical reactions [g] . It fails to converge when bimolecular reactions are 
present |7| . This makes the limit (i) attractive in applications. In Section |U we introduce 
the GCM which converges in the limit (i). 

The paper is organized as follows. In Section [21 we summarize the TRM for coupling of 
structured mesoscopic meshes with microscopic simulations. The methodology for simulation 
of stochastic reaction-diffusion processes on irregular meshes and the implementation of the 
CPM is presented in Section |3l The GCM is introduced in Section |H Using numerical 
examples in Section [5l we compare the computational error associated with the TRM with 
that of the GCM for structured meshes and the CPM with the GCM for unstructured meshes. 
We will then discuss the sources of these errors and ways in which they may be reduced. 

2. The two-regime method (TRM) 

The two-regime method (TRM) [lO|, lU] couples microscopic and mesoscopic subdomains by 



careful selection of the jump rate over the interface from the mesoscopic compartments and 
careful placement of these molecules into the microscopic domain. To date, the TRM has 
been used with mesoscopic subdomains with regular meshes 12|, |9|. The advantage of using 



this technique is that accuracy can be gained in 'regions of interest' Qm without the need 
to run computationally expensive microscopic simulations over the entire domain Q. In this 
section we will briefiy cover the two different simulation paradigms and then discuss how 
these paradigms are combined using the TRM. 

2.1. Microscopic simulation 

The defining characteristic of 'microscopic' simulation techniques for diffusion is that each 
molecule in the system is simulated individually on a continuous domain. In particular, these 
techniques follow the trajectory of each Brownian molecule to a resolution dependent on the 
time steps that are used. For illustrative purposes we consider here a time-driven microscopic 
algorithm. That is, an algorithm with a defined constant time step. Furthermore, we will 
not be considering volume exclusion effects in this manuscript. Each molecule, j, is therefore 
considered to be a point particle at some location in space, Xj(t), at time t. The Brownian 
diffusion of these molecules is modelled by (^. Reactions may take place between these 
diffusing molecules at a particular time step if the reactants are within a given reaction 



radius of each other 24 , 22 



Molecule interactions with boundaries depend on the type of boundary: boundaries can 



be reflective, adsorbing or reactive (partially adsorbing) [6]. Considering that y/DAt is much 
smaller than the local radius of curvature of the boundary, then the boundary is locally flat 
on the scale of relative motion of the molecules in one time step. In the case of absorbing 



boundaries, molecules are removed from the system when they are updated to a position 
outside of the boundary. Since we simulate Brownian motion using a finite time step, we 
have to take into account that a molecule can interact with the boundary during the time 
step [t, t + At] even if its computed position at time t + At is inside the simulation domain. 
The probability, Pm that this molecule-boundary interaction occured within the time interval 
(t, t + At] is dependent on the diffusion constant and the initial and final normal distances 
from the boundary the molecule is found (Axj and Axj respectively) 

P^(Aa;„ Ax^, D, At) = exp (^^^) • (3) 

This probability will also be important when it comes to coupling of microscopic simulations 
with mesoscopic simulations via an interface in the two-regime method 10, 111 



2.2. Mesoscopic simulation 

Mesoscopic approaches to reaction-diffusion processes are simulated on a lattice. For the 
purposes of the TRM we will describe how this is done for a regular cubic lattice. The 
distance between each node is h. In a mesoscopic model, molecules can be thought to exist 
only at lattice nodes rather than existing in continuous space. The state of the simulation 
at any moment of time is defined by a set of numbers describing the copy numbers A/ij- of 
molecules of the i-th type at the j-th lattice point. Considering the diffusion of (non-reacting) 
molecules, the expected state of the system E(A/'jj) is described by the equation: 

^^^ = A 5^ (g.,E(M,,) - g,,fcE(M,)) , (4) 

k 

where Qkj is the propensity per molecule to go from the fc-th compartment to the j— th 
compartment. It is possible to show that for a regular lattice with spacing /i, 

{Di/h"^, if k and j are adjacent lattice points, 
0, if k and j are not adjacent lattice points, 

results in the recovery of the discretized form of the diffusion partial differential equation and 
can therefore be used to approximate a diffusion process on the lattice correct to order h"^. 
The simulation of a mesoscopic react ion- diffusion process usually makes use of event-driven 



algorithms, such as the Gillespie algorithm [17| or the Gibson-Bruck algorithm 16|. We 
shall conceptualize the mesoscopic simulation by considering that when a molecule is at a 
particular lattice point, rather than existing at the node, it is somewhere at random inside 
the compartment belonging to the node defined by the lattice dual mesh |5|. That is, for a 
regular cubic lattice with node spacing h, each molecule which is at a particular lattice point 
is thought to exist inside the cubic compartment of side length h for which the lattice point 
is at the center. It is important to note that the state of the molecule has no specific location 
but rather is thought to exist in a probabilistic sense uniformly over its compartment. 



2.3. Interfacing microscopic and mesoscopic simulations 

Interfacing microscopic and mesoscopic simulations of react ion- diffusion processes using tlie 
TRM has previously been derived for mesoscopic regimes that use regular cubic lattices 



10|,lll[. The TRM is proposed by partitioning the domain ^2 into subdomains ([T]) separted by 
the interface / = d^M^^dQc- In both subdomains, molecules behave as they would normally 
according to the rules of that particular regime. We describe the TRM with an event-driven 
mesoscopic simulation in Qq and a time-driven microscopic simulation with constant time 
step At in Qm- Reactions do not cause any issue within the domain because they occur 
locally. We focus, therefore, on the correct manner in which molecules may migrate over the 
interface /. It is assumed that the TRM is simulated such that y/DAt ~ /i <^ 1. A diagram 
of the numerical TRM scheme using a regular cubic lattice can be seen in two dimensions 



in Figured] A detailed TRM algorithm may be found in the reference 10|]. In order that a 
molecular migration over the interface is smooth with optimally small error, the propensity 
r per molecule to cross the interface / from each adjacent compartment is dependent on the 
parameters h and At. For a regular cubic mesoscopic lattice. 



where D is the diffusion constant of the migrating molecule. The TRM considers that 
microscopic molecules in Qm cease to be microscopic molecules, in principle, when they 
migrate over the interface. Molecules are therefore absorbed by the interface / from Qm 
and placed in the closest compartment in Qq- Equation (^ is used to absorb all molecules 
which interacted with the interface. If this is not used then molecules effectively migrate into 
Qc and back out again without changing from a microscopic molecule to a mesoscopic one. 
This is crucial for coupling of the two regimes as outlined in the derivation in the reference 
10(1 . Furthermore, molecules must be precisely placed in ^m when migrating from f2c. In 



particular, the perpendicular distance x the molecule is placed from the interface into Qm is 
found by sampling from the distribution f{x) 



/(a;) = ./— ^erfcf JL^ ), (7) 

^^ ^ \ 4:DAt VVlDAt;' ^ ^ 

where erfc (x) = S/vr J exp(— t^)(it is the complementary error function. In higher dimen- 
sions, the initial position of molecules migrating into Qm can be chosen to be uniformly 
distributed tangentially to the interface in the region of the originating compartment [1 li] . 
Then the error associated with the TRM is 0{h). We shall investigate the error associated 
with the TRM in ID in a later section of this manuscript and compare it with the GCM 
method introduced in Section HI 

3. Compartment-placement method (CPM) 

In this section, we will discuss how mesoscopic simulation is irnplemented on an irregular 
lattice [5] . We will then present a brief description of the CPM 20 . 
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Figure 1: Graphical representation of the TRM on a regular square lattice. 



3.1. Mesoscopic simulation on unstructured meshes 

Mesoscopic simulations on Cartesian meshes are convenient in the sense that they are mem- 
ory lenient. However, complex geometries and surfaces with high curvature, are easier to 
resolve accurately with an unstructured mesh. Living cells can have different shapes and 
eukaryotes have a complicated internal structure with two-dimensional membranes and a one- 
dimensional cytoskeleton [l| . The geometrical flexibility of unstructured meshes is therefore 
an advantage when considering simulations of realistic biological problems. 

Consider a domain Vt. The domain is covered by a primal mesh, such that the bound- 
ary dVt is covered with non-overlapping triangles and the domain VL is covered with non- 
overlapping tetrahedra (resp. triangles in 2D). A dual mesh is constructed from the primal 
mesh, see Figure 121 from the bisectors of the tetrahedra (resp. triangles) that use the nodes 
as vertices. The diffusion of molecules is now modelled as discrete jumps between the nodes 
of the dual mesh. The rate g^j at which a molecules jump from voxel Vi to Vj is given by 
the diffusion constant of the molecule and the finite element discretization of the Laplacian 
on the primal mesh. For details on how the unstructured meshes and the diffusion matrix 
are generated the reader is referred to |5|. 

3.2. Interfacing microscopic and mesoscopic simulations 

The algorithm for the CPM is presented in a similar way to the TRM. The algorithm pro- 
gresses asynchronously by updates in the mesoscopic simulation and microscopic simulation 



separately [20|]. The jump rates from compartments that are on the interface / between 
regimes are calculated from the underlying mesh over the entire domain. That is, the jump 
rates are calculated by computing the mesoscopic jump rates between interfacial compart- 
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Figure 2: Schematic of CPM computational dom,ain. (a) The prim,al m,esh is indicated with red dashed lines. 
The nodes are connected to form triangles. The bisectors are then drawn in to create the dual mesh (blue 
dotted lines). Compartments are drawn from the dual mesh with one node at the center of each compartment. 
One example compartment is shown in blue. 

(b) The domain is split into mesoscopic Qc o^c^ microscopic Qm domains. Jumps between compartments 
and from the compartments into Qm <ife calculated using a finite element discretization of the Laplacian. 
The copy numbers of molecules in each compartment in Qc o,f^ stored whilst in Qm each molecule has its 
own position in continuous space. 



ments and "compartnients" that are adjacent to the interface in the microscopic domain VLm 
(see Figure [2]). 

Molecules that start in a compartment in VLq and, at the end of the time step, have 
ended up in VLm are initialized uniformly inside the "compartment" which they jump into, 
and is the process from which the CPM has been named. Molecules in VLm migrate back 
to Vtc via microscopic domain diffusion (|2]). When a molecule appears inside one of the 
mesoscopic compartments from Brownian motion, it is encorporated into that compartment 
by increasing the copy number inside this compartment. 

The CPM has been determined using heuristics. Molecules that are in compartments 
obey mesoscopic rules for diffusive migration. This includes molecules that are on interfacial 
compartments. They jump to compartments in VLm as though they were still in Vie- When 
this occurs, initialization of the molecules must take place. The molecules are initiated 
uniformly over the compartment in which they are placed. Molecules are not placed at 
the node at the center of this compartment because this would unphysically concentrate 
molecules at this point and reactions would occur between possible reactants apon migration 
over the interface. Conversely, molecules that are in VLm obey microscopic rules for diffusive 
migration (Brownian trajectory). When this Brownian trajectory leads to a compartment, it 
can no longer be described using the microscopic description and is added to the compartment 



[G.l] Initialize lattice over whole domain Q and construct dual mesh (compartments). 
Generate interface / on the edges of compartments to separate Qm from Qq- Choose 
At and set time t = 0. Determine qkj using finite element method between all 
compartments |5|. 

[G.2] Initialize the initial state of the system by placing molecules in compartments in 
Qc and placing molecules in continuous space in Qm- Count and store numbers 
of molecules in ghost cells, those compartments in Qm which are adjacent to the 
interface I. 

[G.3] Determine the time r for the next event (reaction or diffusive) in Qc or diffusive 
jumps to and from ghost cells and Qc- 

[G.4] If t + r < At + At[t/At\ then change the state of the system to reflect the next 
event corresponding to r and update time t := t + r. If this event is a diffusive jump 
from ghost cell to Qc choose a molecule at random within the relavant ghost cell 
to migrate. If this event is a diffusive jump from Qq to a ghost cell then initialize 
this molecule with uniform probability over the ghost cell. 

[G.5] If t + r > At + At[t/AtJ then update the positions of all molecules in Qm using 



([2]). Check for reactions in Q,m [Hi- All molecules incident on the interface I are 
reflected. Update time t := At + At[t/AtJ 
[G.6] Repeat steps [G.3]-[G.5] until the desired end of the simulation. 



Table 1: The ghost ceh method algorithm. 

in which it lands. As we shall see, this heuristic approach can lead to inaccuracies. The 
inaccuracies can be minimized if h"^ ~ DAt (that is, if the size of the compartment is 
approximately the size of a microscopic molecular jump). 

4. The ghost cell method (GCM) 

Here we will consider a new method for interfacing mesoscopic and microscopic simulations. 
This method uses different assumptions to the TRM and CPM and is therefore implemented 
differently. We call this method the ghost cell method (GCM) since microscopic molecules 
in Qm feel the presence of a pseudo-compartment allowing for instantaneous jumping from 
Qm to Qc in the same way that molecules within compartments jump instantaneously. The 
steps of the GCM are given in Table [1] 

The key assumption that is used in the TRM and CPM is that molecules in Qm mi- 
grate via diffusion (|2]) over the interface / whereby they become parts of the corresponding 
compartment. In the GCM, this assumption is relaxed. Instead, molecules migrate over the 
interface using the compartment-based approach in both directions. Microscopic molecules 
in Qm near the interface feel the presence of a layer of "ghost" cells (compartments). In the 
step [G.2] in Tabled! we calculate the numbers of molecules in these "ghost" cells. They are 
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Figure 3: A diagram illustrating the fundamental differences between (a) the TRM/CPM paradigm and 
(b) the GCM paradigm. 



used in the step [G.4] to create a fully compartment-based simulation of transition across 
the interface /. 

To justify the GCM, let us consider a simulation of diffusion in a domain Q for which a 
mesoscopic method was implemented. Then consider the same domain where a microscopic 
simulation is implemented. Let the molecules of the microscopic simulation be "binned" 
according to compartments of the mesoscopic simulation. The expected number of molecules 
binned into each compartment should match that of the mesoscopic simulation to within 
the precision of the mesoscopic method. This is because both simulations are accurate 
representations of the same phenomena, diffusion. This is the philosophy behind the GCM. 
Molecules which are binned into ghost compartments near the interface may jump into 
compartments in Qq via the rates prescribed by the mesoscopic algorithm. If both regimes 
are correct individually then the flux over the interface / is the same as though a mesoscopic 
algorithm was used over the whole domain. To ensure that microscopic molecules do not 
migrate to fie via diffusion (|2]), they are reflected at the interface / in the step [G.5]. Figure [3] 
demonstrates the principle differences between a TRM/CPM and a GCM description of the 



interface. In Appendix A we provide a mathematical analysis of the GCM in one dimension 
to demonstrate that the expected concentration and flux of molecules over the interface are 
matched. The theoretical error associated with the GCM scales as y/At which is on the 
same order as that of the TRM. Unlike the TRM, this error, as we will see in the later part 
of this manuscript, is reduced to zero by reducing y/At/h. 

The ghost cell method is implemented using the algorithm in Table [1] This algorithm is 
given for an event-driven mesoscopic simulation and a time-driven microscopic simulation, 
however it can also be extended to event-driven microscopic simulations. 



5. Numerical results and discussion 

In this section, we present numerical examples comparing the TRM, CPM and GCM. First, 
we demonstrate how the error associated with the interface I is dependent on choices of 
mesh spacing h in the mesoscopic subdomain at the interface and the time step chosen 
for the microscopic subdomain At for both the TRM and the GCM using one dimensional 
simulations. 

5.1. One dimensional simulations: TRM versus GCM 

We use a simple diffusion test problem to compare the diffusive flow over the interface with 
an exact solution which can be analytically obtained. We use the domain Vt = (0, 1) and 
subdomains flc = (0, 0.5), Qm = (0.5, 1), which are separated by the interface / = {0.5}. We 
initially position A^o = 5 x 10^ molecules according to the distribution g{x) = 2x, x E fl. We 
construct regular spaced compartments of width ho = 0.1 within Qc and "bin" the molecules 
generated in Qc into these compartments. We allow these molecules to diffuse throughout 
the domain Q with a diffusion constant D = 1 using the TRM or GCM until time t = 1. At 
the boundary a; = molecules are absorbed and placed at x = 1. At the boundary x = 1 
molecules are reflected. In this way, Nog^x) is the steady state distribution of this system 
and 0.25A^o is the steady state number of molecules in Qc- We define a measure of the error 
E to this test problem for each simulation scheme 

A^o 

where J\fj is the copy number of molecules in the j-th compartment evaluated at t = 1 and 
the sum is taken over all compartments in Qc- 

In order to see the effect of the compartment spacing near the interface h on the error E 
for both the TRM and GCM we start with the set of regular compartments 

(0,M, iho,2ho), ..., (0.5-/10,0.5), 

which have nodes (centres of compartments) at /io/2, Sho/2, . . . , 0.5 — ho/2. Then we use 
the following lattice refinement technique designed specifically so that the position of the 
interface does not change (see Figure H]) : 

[R.l] Delete the two nodes closest to the interface. 

[R.2] Introduce into the space between the new node closest to the interface and the interface 

(a distance of Aa;) three nodes placed consecutively a distance of 2Ax/7 from the node 

to its left. 
[R.3] Recompute the compartments by finding the bisectors of each node. 

The specific distances in the step [R.2] are chosen such that the interface does not change 
location and the last two compartments have the same size. This is also the size that is given 
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Figure 4: Diagram of one iteration of the lattice refinement [R.1]-[R.3]. 



to the ghost cell in the GCM. The refinement algorithm [R.1]-[R.3] is repeated m times such 
that the size of the final compartment in Qc (and ghost cell), h^ 
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K 



ho 



(9) 



A diagram representing one iteration of the refinement technique [R.1]-[R.3] is shown in 
Figure JH The error is computed for various final compartment sizes hm {m = 0, 1, . . . , 10) 
and various time steps Atk (fc = 0, 1, . . . , 10) where 



Atfc = 2'=Ato, 



(10) 



and Ato = 5 X 10"^ 

Figures [5] and [6] show how the absolute error ||i?|| given by ([8]) depends on both parameters 
hm (compartment size on the interface) and At for the TRM and GCM algorithms respec- 
tively. The error due to the interface in the TRM includes a shift of /irra/2 in the expected 
distribution of molecules at the interface into Qq [10]. This is because of the "initialization" 
of molecules from Q^ into Qc- Unlike the initialization of molecules from Qq into ^m, 
molecules that are transported in the reverse direction cannot be placed carefully according 
to a continuous distribution but must necessarily be placed in the nearest compartment. 
This initialization has an expected position of /im/2 away from the boundary causing a shift 
of /im/2 in the distribution of molecules. However, if molecules could be initiahzed into fie 
with a continuous distribution, for symmetry reasons one would expect this to be done with 
a distribution of f{x) given by (jTj). The average distance, therefore, that a molecule would 
ideally be placed into Qc is /q xf{x)dx = y/7cDAt/2. Therefore, the error that is due to 
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Figure 5: Surface plot of the absolute error \\E\\ defined by ^ as a function of compartment size at the 
interface km and time step At for the one dimensional test problem using the TRM. 



unphysical shifting of molecules is proportional to the expected shift of molecules as they 
are transferred from Hm to ^c- That is E oc hm — vTrDAi. 

In Figure [5] a dotted red line showing hm = vnDAt approximately follows the path of 
the minimum absolute error. The discrepancy between the actual minimum absolute error 
and the dotted red line in Figure [5] can be attributed to higher order error that is inherent in 
the mesoscopic approximation to the diffusion equation. To show that E (x hm — \ArDAt, 
Figure [7] is a plot of error E versus hm — ^A^DAt. The plot is generated by using various 
values of hm (see legend) and then plotting a number of points while changing At. Whilst it 
is clear that the graph is approximately linear, the higher order mesoscopic error is clearly 
seen in the form of a vertical displacement of this curve about the origin. The effect that 
the higher order mesoscopic error has on the interface is difficult to quantify because it will 
depend on the particular molecular system. Therefore, the best choice of parameters that 
can be chosen for the TRM is hm ~ VnDAt. 

In Figure [6], we see that the error of the GCM depends on At and specifically on its 



relative size compared to h (the analysis of the GCM is provided in Appendix A ) . Rapidly 
increasing error (quickly saturating the color bar in Figure[6]) is observed when hm ~ vTcDAt. 
The higher order mesoscopic error artefact can also be seen in Figure M since this artefact 
is independent of the coupling mechanism (see the larger absolute error for large values of 
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Figure 6: Surface plot of the absolute error \\E\\ defined by ^ as a function of compartment size at the 
interface km and time step At for the one dimensional test problem using the GCM. 



h). The GCM is therefore most accurate for very small values of At. Whilst in practice 
making At small may significantly increase the computing time, small At is often required 
for accurate microscopic simulation (for example, capturing reactions with high resolution) 
and in such cases the GCM is more appropriate than the TRM. 

5.2. Three dimensional simulations: CPM versus GCM 

In this section we will demonstrate how, when using an unstructured mesh, the error asso- 
ciated with the GCM coupling converges as At — )■ whereas error associated with the CPM 
is minimized when h ~ y/DAt where h is the average size of boundary compartments. Both 
the error associated with the CPM and GCM are due to imbalances in the flux of molecules 
over the interface. We implement the CPM and GCM in three spatial dimensions using a 
tetrahedral primal mesh as described in Section [31 The implementation builds on the freely 
available software URDME fi]. 

We consider a cube with side length L = 1. The cube is first discretized with an unstruc- 
tured mesh and then divided into a mesoscopic region Qq, and a microscopic region Qm, 
where flu is the set of all voxels with a vertex (x, y, z) such that x < 0.5 and Qq = ^\ ^m- 
Here Q is the set of all voxels. The partitioning is illustrated in Figure M for two different 
mesh sizes. We start each simulation with A^o = 2 ■ 10^ molecules whose initial positions are 
sampled from a uniform distribution. The diffusion constant of the molecules is D = 1, and 
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Figure 7: Scatter plot of the error E versus ft,,„ — y/wDAt for the TRM. The different color points represent 
different values of the compartment size at the interface h (see legend) and in each instance At is varied 
from 5 X 10"^ to b x 10~3_ 



we simulate the system for time t=0.1. Since we start with a uniform distribution and the 
molecules only diffuse and do not react, we expect the distribution to be uniform at the 
final time. As the interface is parallel with the y — 2;-plane, we expect that the distributions 
of molecules in the y- and 2;-directions are uniform, but that we get a small error in the 
distribution of molecules in the x-direction. We now divide the x-axis into 10 bins of equal 
length, and then count the number of molecules in each bin at the final time. Mesoscopic 
molecules are binned by first sampling a continuous position from a uniform distribution on 
the voxel. We expect Nq/IQ molecules in each bin, and can therefore estimate the error E 
by 



E 



E"i \N^ - iVo/10| 



N. 



;ii) 



In Figure |9] we have computed ||ii^|| for different mesh sizes and time steps. As expected, the 
error decreases as we refine the mesh and decrease the time step. 

In the CPM method, mesoscopic (resp. microscopic) molecules stay mesoscopic (resp. 
microscopic) during a time step. This implies that the time step should be chosen sufficiently 
small such that a molecule does not diffuse across several voxels. On the other hand, if the 
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Figure 8: The cube [0, 1]^ is partitioned into a mesoscopic region (grey) and a microscopic region (white). 
(a) a coarse mesh; (b) a fine mesh. 



10" 



■25393 voxels 
■49101 voxels 
■72413 voxels 



LU 




Figure 9: The error \\E\\ of the GCM method for different mesh sizes and time steps. The error decreases 
with decreasing time step, as expected. 
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Figure 10: Comparison of the error \\E\\ produced by the GCM and CPM for interfacing microscopic and 
mesoscopic simulations as a function of the time step in the microscopic simulation domain Ai. The length 
scale h is defined to he the cubic root of the average volume of a voxel. 



time step is too small the distribution of molecules in space will be biased towards the 
microscopic region. This can be seen by considering a microscopic molecule diffusing into 
the mesoscopic regime. If the time step is small, it is likely that it will be close to the 
microscopic regime at the end of the time step, but if it ends up on the mesoscopic side it 
will nevertheless be considered uniformly distributed in the voxel at the end of the time step. 
Thus, the time step should not be chosen too small relative to the size of the voxels, or the 
error due to the spatial splitting will become large. 

Since the GCM converges with decreasing time step, but performs worse for larger time 
steps, one could suspect that there is a regime where the CPM in 20|] performs better than 
the GCM. At some point, however, the error of the GCM will become small and outperform 
the CPM in |20| . The errors of the different methods are compared in Figure [10] for a mesh 
with 49101 voxels. Indeed, we see that the CPM method performs better for time steps 
down to almost At = 10~^, at which point the error of the GCM method becomes smaller. 



6. Summary 

In this paper we have compared two existing mesoscopic-microscopic coupling techniques for 
stochastic simulations of reaction-diffusion processes with a new convergent method called 
the ghost cell method (GCM). Here we will summarize the specific sources of error of the 
TRM, CPM and GCM, when they converge, how they may be optimized for accuracy and 
notes on their computational efficiency. 
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6.1. Summary of the two-regime method 

The TRM couples molecules by considering that mesoscopic compartments contain molecules 
that are evenly distributed in a probabilistic sense. As these molecules diffuse over the in- 
terface they are placed according to the distribution f{x) given by ([7]). Molecules migrating 
in reverse from the microscopic regime to the mesoscopic regime must be absorbed by the 
interfacial compartments and be indistinguishable from other molecules in these compart- 
ments by virtue of this paradigm. As such, instead of migrating an average distance over the 
interface proportional to \/At it becomes evenly distributed over the compartment with an 
expected location of h/2 over the interface. The molecules are therefore effectively shifted 
{h — \/TxD/S.t)/2 into the compartment regime. This shift in the molecules therefore creates 
a discontinuity of in the distribution to find molecules on the interface and therefore an error 
due to the presence of the coupling proportional to h — yTiDKt. The nature of this error is 
that, if the expected net flux of molecules over the interface is 0, then no error due to the 
presence of the interface will be experienced. This is important to note, since, this is not the 
case for both the CPM and GCM methods. Furthermore, since the error is proportional to 
h — y/iiDAt it clearly converges in the limiting case (ii) described in the introduction but 
not in the limiting case (i). 

Whilst the TRM can give controllably accurate results, it can be computationally more 
costly to implement. This is because perfect absorption of molecules is required on the 
boundary and this means that each molecule in the molecular-based domain needs to be 
checked for interaction with the boundary in a given time step using (E]). 

6.2. Summary of the compartment-placement method 

The CPM is a coupling mechanism that, whilst heuristically derived, can produce accurate 
results under some circumstances and do so with minimum computational cost. Molecules 
are placed within a pseudo-compartment in the molecular-based domain via diffusion from 
the compartment-based domain. In reverse molecules are placed in compartments from 
the molecular-based domain via diffusion of these molecules in the continuous domain over 
the interface. The antisymmetry that is seen in the methods of migration, mesoscopic 
to microscopic and microscopic to mesoscopic, result in a boundary layer in the expected 
distribution of molecules at the interface. This boundary layer is caused by the fact that 
molecules diffusing from the molecular region to the compartmental region are considered 
uniformly distributed on the compartment at the end of the time step. If At is small 
compared to h, this will be a poor approximation. It should thus be noted that the error of 
the CPM does not converge in the limiting case (i) described in the introduction, however 
the error appears to converge according to limiting case (ii). 

The CPM is computationally efficient. Its only inefficiency is that, unlike the TRM, it 
requires the knowledge of a pseudo-compartment in the microscopic domain. The TRM is 
therefore more appropriate than the CPM for coupling completely independent simulation 
algorithms, since the CPM requires its own custom algorithm to be implemented fully. 
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6. 3. Summary of the ghost cell m,ethod 

The GCM couples molecules according to a discrete domain on each side of the interface. 
Molecules that are in the microscopic domain are binned according to a ghost compart- 
ment/cell and jump into the mesoscopic domain using jump rates derived using the meso- 
scopic approach. In such a way, symmetry is conserved in the method of migration from 
mesoscopic to microscopic and from microscopic to mesoscopic, unlike the CPM. It is im- 
portant that the molecules are binned correctly for this coupling to work accurately. The 
error, therefore, can be attributed to molecules that are in the ghost cell when they should 
not be, or not in the ghost cell when they should be. Therefore, if the compartment size 
h at the interface (and of the ghost cell) is much larger than the resolution of the particle 
tracking in the microscopic domain, the correct number of molecules will be in the ghost 
cell. The error therefore converges in the limit of small At so long as h is sufficiently coarse. 
Furthermore, it is possible to show that, unlike the TRM, this source of error is not due to a 



displacement of molecules but an unballanced flux of molecules (see Appendix A ) and will 
therefore appear even if the expected net flux over the interface is 0. The GCM, however, 
is convergent in the limiting case (i) but not (ii) from the introduction giving the GCM a 
unique advantage over both the CPM and TRM. 

The GCM is computationally efficient for small At since the jump rates from the micro- 
scopic domain to the mesoscopic domain are determined by the ghost cell size and not the 
time step (like the TRM for example). 
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Appendix A. Mathematical justification for the ghost cell method 

Here we present an analysis of the GCM in one- dimension. We show that error of the 
GCM that is produced on the interface between mesoscopic and microscopic subdomains 
converges in the case (i). Specifically, we see convergence of the interface-derived error as 
A = yDAt/h — )■ 0. This property of convergence is unique to the GCM when compared with 
other reported coupling mechanisms. In showing that the interface-derived error vanishes in 
the small time step limit, we will show that rapid variation within the boundary layer of the 
interface vanishes as A — )■ 0, leaving behind a linear approximation of the true distribution 
of molecules. Since the error at the interface will be of order h'^ it is accurate to the same 
order as the mesoscopic algorithm itself. 



Without loss of generality, consider an interface at x = on an infinite one-dimensional 
domain. To the left of this interface (x < 0) a compartment-based model is used with fixed 
compartment size h. To the right of the interface (x > 0) a molecular-based algorithm is 
used and is updated at fixed time increments of At, i.e. flc = (— cxd,0) and Qm = (0,cxd). 
We denote the compartments in Qc by Ck = {—kh, —kh + h), where k = 1,2, . . . . Then the 
interface compartment is Ci = {—h, 0). The ghost cell will be denoted by Cm = (0, h). 

Molecules in flc are described only by their compartment. Their compartment changes 
with an exponentially distributed random time with a rate that is given by D/h"^. This rate 
is conditional on initial and final states being compartments. The rate given by D/h^ is 
chosen in such a way that the expected number of molecules in each compartment matches 
that of a discretized diffusion equation (see (jl]) and (Q). These rates, however, breakdown 
in the case of the TRM because the initial and final states of jump across the interface 
are not compartments but rather the final state is a molecule in Qm- The jump across the 
boundary for the TRM is given by iQ. Molecules in Qm have one thing in common with 
those in compartments. A domain that is modeled microscopically and then binned into 
compartments shows the same expected behaviour as a domain modeled with compartments 
to leading and first order accuracy in h in the limit as At —t- 0. Therefore, in an attempt to 
interface the two regimes together it may be appropriate to bin molecules in Qm into a ghost 
cell/compartment Cm near the interface. The molecules that are in Cm will then have the 
same properties as a compartment from the perspective of the interface compartment Ci. To 
this end, any molecule in Cm may spontaneously change state from the molecular domain 
to Ci. We expect that since the interface compartment-bound molecules see a compartment 
state for molecules in Cm, the change of state from Ci to a random position within Cm will 
occur with a normal inter-compartmental rate. In the following analysis we show that this 
is the case. 

We shall test the hypothesis by matching the master equations for Ci and the probability 
distribution in Qm in such a way that no rapid variation in probability to find molecules, 
p{x,t), is apparent at the interface. We shall assume that the rate for molecules to jump 
into Qm from Ci is F"*" and are placed in an initial position from the interface given by 
the probability distribution f{x). Molecules in Qm spontaneously jump into Ci with a rate 
r~g{x). Functions f{x) and g{x) are normalized such that they have a unit integral over 
Qm- We shall show that 

r^ = r- = ^ and 9{x) = f{x)=h fo^^^^*^' (A.l) 

I 0, otherwise. 

We will find it convenient for the sake of notation to introduce the parameters 



, h^T^ , , /DAt 

a = — — — , and A = — ; . 

D h 



To show ( 1A.1I) we focus on the purely diffusive problem, since bulk reactions have no effect on 



boundary conditions. In order to limit the flux of molecules jumping into Ci, all molecules 
in Qm that hit the interface by Brownian motion are reflected back to Qm- 
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We denote the probability of finding a molecule in compartment Ck, k = 1,2, ... , by 
Pk{t)h (so that Pk(t) approximates the probability density function at the node within this 
compartment). If we denote by p{x,t) the probability density function of the discrete- 
time molecular-based algorithm, then the transmission/refiection rules give us the following 
master equation 



roc 

Pi{t + At) = (l-(l + a+)A2)pi(t)+A2p2(t) + a-AM g 

Jo 



{x)p{x,t)dx, (A. 2) 



p{x,t + At) 



p{y,t) 







x-y)\ , f-{x + y) 



^^P — fta — + ^^P 



y/AnDAt [ V 4DAt J ^ \ AD At 



dy 



a+DAtp^{t)f{x) g- DAtp{x,t)g{x) 
+ h h • ^^-^^ 



In the vicinity of x = there is a boundary layer of width 0{h) so long as f{x) and 
g{x) vanish for x ^ /i [6|. We rescale (]A.2p and (JA.31) using the (dimensionless) boundary 
layer coordinate ^ = x/h. We also denote Pinner{iii) = p{h^,t), fmneriO = hf{h^) and 
fi'innerlO = hg{h^). The rescalings of / and g hj h are done to keep the integrals of these 
functions equal to 1. Thus, in the boundary layer coordinates, (1A.2I) and ( 1A.3I) become 

Piit + At) = {l-il + a^)A^)pi{t)+A^p2{t) 

+ a- A'' / (7inncr(0 Pinner 1^, t) d^, (A.4) 

Jo 

/•oo 

Pi„ner(e, t + At) = A'' / pi^,,,,{r], t) [K {A'^t] - O) + ^ {^'\v + O)] dr] 

+ a^A'^Pi(t) /inncr(0 " a~A^Pinner(^, t) ginner{Oy (^-5) 

where i^(x) = (47r)~^/^ exp(— x^/4). The parameter A gives us the relative size of DAt to 
h'^ and we wish to show that as A — )■ the distribution of molecules across the boundary is 
smooth and the error that remains is of the order of h"^, which is the same size of the error 
associated with the mesoscopic discretization in Qc- In order to join these models smoothly 
we require in Qq that 

Pi(t) = p(-h/2,t)=p{0,t)-^p,{0,t) + O{h') + O{A), (A.6) 

P2{t) = pi-3h/2,t)=p{0,t)-—p,{0,t) + O{h') + O{A), (A.7) 

Piit + At) = pi(t) + 0(At), (A.8) 

while, for the molecular-based side, we want variation from the linear approximation in the 
boundary layer to be hmited to 0(A) up to order h"^ accuracy, so that 

Pinneri^,t) = ^(0, t) + /l ePx(0, t) + ©(/l^) + 0(A), (A.9) 

Pinner{^,t + At) = Pi^,,,, {^ , t) + O (At) . (A.IO) 
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The prescription of a consistent probability density p{0,t) and derivative Px{0,t) for both 
sides of the interface, along with linear approximations sufficiently close to the interface 
equates to continuity and differentiability over the interface which is the matching condition 
that we are attempting to achieve. 

Substituting flA.6p - flA.10|) into flA.4|) and flA.Sp and equating terms of the same order in 



h and leading order in A gives the following conditions that must be placed on g{x), f{x), 
a~^ and a~: 

(i) 0{h^A^) terms from equation flA.4p give condition: 

POO 

a+ = a- / dinner (0,de (A. 11) 



Condition (IA.11|) states how the relative rates for molecules to transition to and from 
Qm and Qc niust be dependent on the relative sizes of Ci and Cm- 



(ii) 0{h^A^) terms from equation flA.4|) give condition 



2 = a+ + a- e^inner(0,de- (A.12) 

Condition ( JA.121) states how the rates for molecules to transition to and from Qm and 
Qc depend on the average distance molecules are placed from the interface when placed 
within Cm- This is the same condition given in Qc for jumps between the compartments. 

(iii) 0(/i°A°) terms from equation (IA.5|) give condition: 



a^/mner(0 = " fi'inner(0- (A. 13) 

Condition (IA.13|) states that molecules must be placed into Qm with the same probability 
weighting that they are taken out and placed back into Qq- 

(iv) 0{h^A^) terms from equation (lA.SP are automatically satisfied. 

The GCM that is presented in this manuscript uses parameters (lA.ip which satisfy the three 
conditions flA.ll|) - flA.13|) listed above. Such a scheme, therefore, has an error that is no 
greater than the error of the mesoscopic scheme in the limit A — )■ 0, in other words, in the 
limit At — )■ whilst h remains constant. 
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